Maximum velocity of a fluxon in a stack of coupled Josephson junctions 
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Dynamics of a fluxon in a stack of inductively coupled long Josephson junctions is studied 
analytically and numerically. We demonstrate that the fluxon has a maximum velocity, which 
does not necessarily coincide with any of the characteristic Josephson plasma wave velocities. The 
maximum fluxon velocity is found by means of numerical simulations of the quasi-infinite system. 
Using the variational approximation, we propose a simple analytical formula for the dependence of 
the fluxon's maximum velocity on the coupling constant and on the distribution of critical currents 
in different layers. This analysis yields rather precise results in the limit of small dissipation. The 
simulations also show that nonzero dissipation additionally stabilizes the fluxon. 

74.50.+r, 41.60.Bq, 74.80.Dm 



I. INTRODUCTION 



Experimental and theoretical studies of magnetic flux quanta (fiuxonsl m stacks of inductively coupled long Joseph- 
son junctions (LJJ's) have recently attracted a great deal of attentiorJjQ. The interest to the stacks is stimulated 
both by the fact that the high-T c superconductors, on the atomic level, have a naturally layered structure that is 
tantamount to an intrinsic Josephson stacktl, and by the development of the (Nb-Al-A10 x )jv-Nb low-T c technology . . 
which is demonstrated fabrication of artificial stacks of up to 28 LJJ's, with a parameter spread between them < 10%L. 

A fundamental characteristic of a single LJJ is its Swihart velocity cq, i.e. , a minimum phase velocity of the 
electromagnetic (Josephson plasma) waves propagating in the superconducting micro-strip structured. Simultaneously, 
Co is the maximum velocity for fluxons that correspond to topological solitons of the sine-Gordon model describing 
LJJ. In a system of N linearly coupled junctions, the dispersion curve has N branches corresponding to different 
modes of the linear electromagnetic waves propagating in the system. Accordingly, there are N split (different) 
Swihart velocities n = l,2,...,N (see, Refs. | 10) such that c„ < c^S. For example, in the simplest case of 

two coupled junctions, there are two Swihart velocities c\ = c_ and c 2 = c+ (c_ < c$ < c+), which correspond, 
respectively, to the system's in-phase and out-of-phase Josephson plasma wave eigenmodes. In the AT- fold stack, we 
also use notation c_ = c£ and c+ = which are the smallest and largest Swihart velocities. 

An important issue is to study the conditions for the existence and stability of a single-fluxon state in the stacked 
system. It is implied that the fluxon is trapped in one junction and its screening currents spread over neighboring 
junctions. The fluxon induces, through the magnetic coupling, "images" in adjacent layers, so that a full solution for 
the single fluxon state in the stack includes both the core topological soliton in the central layer and its non-topological 
images in the other layers. Such a fluxon state we denote as [0| . . . |0|1|0| . . . |0]. 

As mentioned above, in a standard singleJaarrier LJJ a fully stable fluxon with a velocity exceeding th&nSwihart 
velocity cannot exist. It was first suggestedtil, and recently demonstrated theoretically and experimentallyalla, that 
a fluxon may nevertheless move in a multilayer system with a velocity which exceeds the minimum phase velocity of 
the plasma waves. It is important to find the maximum fluxon velocity u max and its dependence on the parameters 
of the system. A possibility of having M max > c_ is especially interesting, as it implies steady motion of the fluxon at 
u > c_ with Cherenkov radiation tail of Josephson plasma waves behind it. This issue, which is of evident physical 
interest, is the main subject of the present work. The possible range of u max was not investigated systematically in 
Refs. ^,^[ The only prediction which has been made is that u max > c_ for asymmetric junctions in a two-fold stack 
(for the case when the fluxon is trapped in LJJ with lower j c ), and that u max > c_ always holds in A^-fold stacks of 
identical junctions. The exact value of it max has not been found until now. 

We will discuss three cases which differ by the number N of the coupled LJJ's. These cases are N = 2, 3, and 
oo. For N = 2 and N = 3 we will consider a system of asymmetric LJJ's with different critical currents j c , in which 
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the maximum velocity is different depending on where the fluxon is trapped. In the case N — 3 and N — oo we will 
assume that the core topological soliton is placed in the central junction which will be labeled by 0, while LJJ's above 
and below the central one will be labeled by 1,2, . . . , oo and —1, —2, . . . , — oo, respectively. 

In section || we formulate the model, section III displays results of full PDE numerical simulations of the asymmetric 
model for the cases N — 2 and N = 3. To choose an analytical form for the fitting function which predicts the 
dependence zt max on junction parameters, in section IV we use the variational approximation (VA). Although VA 
does not produce very accurate quantitative results, it predicts reasonable functional dependence for u max . Section 
^ concludes the work and summarizes the obtained results for different N. 



II. THE MODEL 



A model for iV-fold stack of long Josephson junctions is well know 

(0n) xa; = [<t>n)tt + sin (p n + 7 + a((f>n)t 

-S[(<f> n -i) tt + sin </>„_! + a(<p n -i)t + (<Pn+i)tt + sin0„+i + a(<p n+ i) t + 2j] , (1) 

where <p n is the Josephson phase across the n-th LJJ, n = —N/2 . . .N/2, the coordinate x and time t are measured 
in units of the Josephson length Aj and inverse plasma frequency lu^ 1 of single-layer LJJ, 5 < is a dimensionless 
coupling parametenij, a is a dissipative constant, and 7 is the density of the bias current flowing through the stack. 
We consider the most natural case when the bias current is the same in all layers. 

The model (|l|) pertains to a stack consisting of an infinite number of junctions, or to an exotic configuration, in 
which the stack is closed in a loop in z direction (0at+i = (pi). In practice, the equations for the edge (top and 
bottom) junctions include only a half of the coupling terms corresponding to the neighboring LJJ's. Note also that 
the model implies that all the junctions are identical. In particular, the Swihart velocity of each uncoupled LJJ, is 
So = 1 in the notation adopted hereafter. 

In the case of N = 2, a releyvant version of the model (|l|), which takes into account the difference in the critical 
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currents of the LJJ's, writes as 



1-5 2 ^°' u ""^ u 1-S 2 
sin cpi S((p ) 



sin 0o- , =«(fo)t+7i (2) 



i)«-^--r^r = «(<Mt + 7, (3) 



1-S 2 ^ Ht J 1-S 2 

where J — j c o/jd is the ratio of the critical currents of the two junctions. When considering the model @ and (|| ) 
below, we will place the fluxon into the LJJ whose phase is denoted as <po. 

Discussing the case of N = 3, we impose the symmetry condition (pi = <p—i, which is natural when the fluxon moves 
in the middle layer. Thus, we can write Eqs. (Q) in the form 

1 _ 1 2|g2 (Mxx - (Mtt - - 1 _ 5 2|g2 (Mxx = 7- a {(pi) t ; (4) 

1 25 
!_2g2 ('Mxx - (M tt - sin0 o - 1 _ 2g2 (cpi) xx = 7- a (0 o ) t ■ (5) 

Note the factor of 2 in the last term of the left-hand side of Eq. (^). 

For further theoretical treatment of the case N = 00, we will assume that the coupling parameter 5 is small. This 
assumption will allow us to write a Lagrangian corresponding to the dynamical equations (|lj) , which is a key ingredient 
of VA. In the case of small 5, neglecting terms ~ S 2 and smaller, one can easily transform Eqs. (Q) into a simplified 
form: 

(<Pn) xx ~ ((pn) tt -sin(P n - S[{(j) n - 1 ) xx + ((p n+1 ) xx \ =a((p n )t-l , (6) 

which we will use for further analysis of the N — 00 case. 

A fluxon steadily moving at a constant velocity u can be described by a solution to the above equations (without the 
a and 7 terms) which depends on the single variable £ = C(x — ut). The constant C is introduced for renormalization 
purposes and will be different in the cases N — 2, 3, and 00. Substituting 



< ^/i-^Cr-ut) (7) 
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into Eqs. (^) and (^|) and neglecting the a and 7 terms (which must be kept if one aims to find an equilibrium velocity 
determined by the balance between the losses and bias current, that is not our objective in this work), we get: 

ff (a) $> +<t>" -sin O = , (8) 
a y >(p 1 +0 O — = . (9) 

The parameter (the subscript 2 implies that the definition is adjusted to the case N = 2), that will be used 
instead of the velocity, is defined as: 



(2) = i-(i-s 2 h 
-s 

For the 3-fold stack, introduction of the traveling coordinate 



a m s - v- . (10) 



transforms Eqs. (W) and (feh into 



r 1 — 9 

5-(«-«t) (11) 



a( 3 V'/ + 0o'- = = O , (12) 
ct (3) 0o + 20i -sin 0o = O , (13) 



where 



r(3) 



_ 1 - (1 - 2S 2 ) u 2 



(14) 



-S 

And, finally, for the case N = 00 the substitution of 

£=(x-vt)/V=S (15) 



into Eqs. (O) yields 



where 



a^Wn + €-1 + €+1 - sin0„ = 0, (16) 



roc) = !-" 2 
-S 1 



(17) 



Thus, from the mathematical viewpoint, the issue is to look for solutions of Eqs. (||) and (|j|), or ( |l2] ) and (|l3|), or 
of Eqs. ( |l6| ) that describe the stationary fluxon. The eventual objective is finding the maximum velocity u max , or the 
minimum value of the parameter a, beyond which the fluxon solution does not exist. The fluxon solution is defined 
by the following boundary conditions: 0o(— 00) = 0, 0o(+oo) = 2ir, and n ^o(±°°) = 0. 

To conclude this section, we recall that the set of the split Swihart velocities can be found in an exact form 
the linearized version of Eqs. (Q), setting there 7 = a = 0. These velocities are given by the following expressior 

cj N) = 1 n = l,2,...,JV. (18) 

/l_ 25cos ( ™ 
V \N + 1 

It is also useful to find the values of the parameter a corresponding to the minimum velocity c_ and the maximum 
velocity among all the velocities given by Eq. @J|). From Eqs. ©, (jHJ) and @, using Eq. dl8|), we find 

(19) 
(20) 
(21) 





= Tl 




= T^2 


a (oo) (c±) 


= T2 
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while zero velocity corresponds to cr 1 - 2,3 ' 00 - 1 = —1/5. 

We can immediately find the values of a m i n for some specially selected values of J. Let us set — 1 in Eqs. (||) 
and (0). This reduces Eqs. (||) and (0) to a simple algebraic equation: 



sin (f> 



sin 0i 
J ' 



(22) 



The state with a fluxon only in the first junction (</>i) can be realized only for J < lilll. This means tha10 

<^(J=1) = 1. (23) 
In the same fashion, setting = \/2 in Eqs. jl^ ) and (|l3|), we get 
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This means that 



sin <po = —j- sin (pi . 



a£t(J = V2) = V2. 



(24) 



(25) 



A similar relation can be obtained for any TV-fold stack. Using Eqs. (||), we get the following result. The value 

la: 

\ 



a min ^ m the stack consisting of (2N — 1) LJJ's is given by the maximum eigenvalue of the N x N matrix, 



( ° 1 
' -1 0-1 



V 



-1 -1 
-2 / 



(26) 



This matrix has the size N x N due to the symmetry which, as we suppose, is present when the fluxon moves in the 
middle layer. In this case 2N — 1 coupled equations can be reduced to N equations in the same way as we reduced 
3 equations for N — 3 to only 2 independent equations. As a result of this reduction, the N-th equation, which 
describes the junction containing the fluxon and corresponds to the last row of matrix contains the factor of 2. 

The values of er m m calculated for some N are summarized in Table Q. One can see that as N — > oo, cr m in — * 2. 
Similar to the above considerations, the single-fluxon state exists only if J < cr m in- Note that this is, actually, quite 
a noteworthy result, which means that, for any N > 2, the steady motion of the fluxon accompanied by the emission 
of the Cherenkov radiation is possible at J — 1, i.e., in the uniform stack of identical LJJ's. 

Another simple case is obtained by setting J = 0. In this case Eq. (^) or Eq. ( |l2| ) yield <fii = 0, hence, the remaining 
equation [(g) or (|l3|), respectively] is nothing else but the single sine-Gordon equation, which has proper solitonic 
solution (giving the [1|0] state) only for positive a. Therefore 



x(0) = 0, for any N . 



(27) 



Physically this means that at large disbalance of critical currents the maximum velocity of fluxon is approaching 

1 EE CQ. 



III. NUMERICAL RESULTS 



In the region u > c_, i.e. at fluxon velocities larger than the lowest phase velocity of plasma waves, the phase 
dynamics in a stack is quite complex. Analytically, it can not be reduced just to looking for solutions of unperturbed 
equation in the form cj> = <p{x — ut). Therefore, direct numerical simulations are necessary to get an insight into the 
problem. 

The PDE's (§) and @ for N = 2 or (|) and (g) for TV = 3 were solved numerically using an explicit method 
[expressing <f> A - ' (t + At) as a function of <f> A,B (t) and (f> A,B (t — At)], and— treating <^> xx with a five-point, <ptt with 
a three-point and (f>t with a two-point symmetric finite-difference schemeEJ'lU. Equations were supplemented by the 
periodic boundary conditions, 4>o{x + t) = <t>o(x) + 2ir, and (j>i(x + 1) = <t>i(x), with the period I = 200, so that 
the fluxon moves in the quasi-infinite system. Numerical stability was checked by doubling the spatial and temporal 
discretization steps Ax and At and checking its influence on the fluxon profiles and current- voltage curves (IVC). The 
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values used for the simulations were Ax = 0.025 and At = 0.00625. To calculate the voltage in each point of the IVC, 
the instant voltage was averaged over the progressively increasing time intervals and also averaged over the length 
of the system. The following convergence criterion was adopted: the difference between the voltages obtained by 
averaging over two successive time intervals must not exceed SV = 5 • 10~ 4 . When the average voltage corresponding 
to the current 7 is found, the current is increased by a small amount £7 = 0.001 to calculate the voltages at the next 
point of the IVC. We use the phases (and their derivatives) attained in the previous point of the IVC as the initial 
conditions for the next point. By gradually increasing 7 and, thus, u, we encountered a maximum value itmax of u 
above which the single-fluxon mode becomes unstable, and the system switches into a resistive state. 

The simulations of IVC's were carried out for the damping values a = 0.02, 0.04, 0.1, in order to understand the 
effect of the dissipation on w max - In experiment, at low temperatures, a can be about 0.01 and less. It turns that 
the simulations of the system with a < 0.02 incurs unaffordablc computational expenses. Therefore, we focus on 
higher values of a and then will try to extrapolate u max (if possible) to the ideal case a = 0, which is of fundamental 
theoretical importance. 

Here, we will only display the numerical results obtained for the cases N = 2 and N = 3. Examples of IVC's for 
N = 2 and 3 and different values of J are shown in Fig. [l] and Fig. |[ Strong effect of the ratio J between the critical 
currents of the junctions on w max can be learned from Fig. ^. In the case J = 2, i.e. , when the fluxon is located in the 
junction with the larger critical current, the fluxon's maximum velocity is smaller than the lowest Swihart velocity 
c_. In addition, one can observe a back bending of the IVC close to its tip. The last point on back bent IVC marked 
as [1| + 1, — 1] corresponds to the state when fluxon-antifluxon pair stretched out in the idle junction with lower j c . In 
the case J = 0.5, i.e., when the fluxon moves in the junction with the smaller critical current, the fluxon's maximum 
velocity is larger than c_. We note that the latter case is nontrivial and of the particular interest since, as it was 
already mentioned above, the fluxon can propagate stably, emitting the Cherenkov radiation. 

Fig. H shows IVC's of the [0|1|0] state for J = 1 and different values of a. As long as we consider 3-fold stack with 
J < v*j the IVC bends to the right into the Cherenkov region. The fluxon motion in such a case is very similar to 
the case J < 1, N = 2. Fig. || demonstrates that the damping not only changes the slope of the IVC at low velocity, 
but also affects u max . 

The numerically found values u max for different values of J are transformed into the corresponding values of the 
parameter a m i n according to Eqs. ( |lo|) and (|l4|). The values of cr m j n for N = 2 and N = 3 are plotted in Fig. || 
and Fig. ^, respectively. Note, that the numerically obtained dependence <7min( J) is in agreement with our analytical 
predictions given by Eqs. (23), (25) and (|27|). These three conditions work better for small a which is quite natural 
since they were derived from the unperturbed (a = 7 = 0) equations. We remind that predictions (^), (|2^) and ( pTj ) 
are strict results and must be considered exact in the framework of inductive- coupling model. We found it interesting 
that the dependence of the maximum current 7 max = 7(u ma x) on J is almost linear with 7 max /J = 0.476, 0.531, 0.623 
for a = 0.02, 0.04, 0.10, respectively. 

In the following section we suggest a functional dependence for cr m i n (u) and will compare this dependence with 
numerical results presented in Fig. || and Fig. ||. 



IV. THE VARIATIONAL APPROXIMATION 



The purpose of this section is to find an analytical dependence which describes the simulation data obtained in the 
previous section. We will present an analysis using VA only for N = 2, but other cases can be considered similarly. 

VA is based on the fundamental fact that Eqs. (p|) and (Q) admit the variational representation with the Lagrangian 
L = jO-2 <i£, where the Lagrangian density corresponding to Eqs. (^]) and (|9|) is 

^2 = ia(^) 2 + ia(0i) 2 + 0^i + (l-cos0 o ) + j(l-cos0!) . (28) 

A crucial step in the application of VA is to adopt an ansatz, i.e., a trial form of the solution. The ansatz is then to be 
inserted into the corresponding Lagrangian, and the integration over £ given by Eq. (Q) must be performed explicitly. 
This will produce an effective Lagrangian as a function of free parameters that the ansatz contains. Finally, the values 
of this parameters will be found from the condition that they must realize an extremum of the effective Lagrangian. 

Thus, it is necessary to select a tractable ansatz (that admits an analytical calculation of the integrals in the 
expression for the Lagrangian), which satisfies the above boundary conditions for the fluxon. Here, both when 
selecting the ansatz and considering the boundary conditions, we ignore the above-mentioned nonvanishing Cherenkov 
oscillatory tails (note that the tail is formally infinitely long in the ideal model with a = 7 = 0, while in the damped 
system the tail decays exponentially far from the fluxon's "body"). The simplest and, in fact, the only practical 
ansatz is 
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cj) (£) = 4arctanexp(A0 , (29) 

^ = B^ML . (30) 
* cosh 2 (AO 

Here, A and B are free real parameters. Note that A can always be defined positive, while B may be positive as well 
as negative. The expression (^0|) was chosen in such a way that it describes the "image" profile qualitatively well in 
comparison with simulation and analytical resultald. 

Still, the effective Lagrangian cannot be immediately calculated with the above ansatz. A further necessary simpli- 
fication is to assume that the amplitude B of the fluxon's images in the noncore layers is sufficiently small, so that the 
nonlinear term sin</>! in Eq. (M) may be linearized. In other words, the term (f — cos^!) in the Lagrangian density 
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(p8|) is to be replaced by \4>\ 



The final expression for the effective Lagrangian is 



l -^ x + J + ^ b2x + Itx + I bx - (31) 

The variation of (|Tj) in B leads to an equation that allows us to eliminate B, i.e., 

fOJA 2 

B (5 + 7Ja\ 2 ) ■ (32) 

The remaining algebraic equation for the fluxon's inverse width A is 

7J 2 a (21a 2 - 5) A 6 - 3 J [25 - 7a 2 (f - 7 J)} A 4 - 15a (14J — 5) A 2 — 75 = . (33) 

This equation has one or three positive roots for A 2 . The first two solutions exist and disappear simultaneously at 
some value of a. The third solution disappears diverging at cr m i n = \Jhj21. This can be seen from the form of the 
coefficient in front of the term ~ A 6 . We drop the third solution from the consideration since it gives a constant 
c m in(>/) which is unphysical due to Eqs. ( p3| ) and (|27|). 

To find the minimum value of a at which two smallest positive roots disappear, we write two conditions: the 
function /(A) ([33]) touches A-axis and its derivative /'(A) vanishes. These two equations determine er m i n and A(<7 m i n ) 
for given value of J. Eliminating, consecutively, the terms ~ A 6 , A 4 and A 2 yields an equation which determines the 
dependence of cr m i n on J: 

56 (7 J + 5) 3 a^ in - (3675 J 2 + 36750 J + 1875) a^ in + 7500 J = , (34) 

A solution to this equation is 



2 36750J + 1875 + 3675 J 2 ± f 25\/l5y (15 + 7 J) (f - 7jf 
CTmin = U2(5 + 7J) 3 ■ (35 ^ 

One can see that this solution exists only for J < 1/7. For J > f/7, Eqs. ( |33|) has only one (third) root, which is 
unphysical. 

To relax the latter limitation, we can make use of the fact that, as it was said above, the asymmetry parameter J 
always takes values J < 1. We will make use of this, treating J as a small parameter, which is also justified by the 
fact that the stack with a sufficiently strong asymmetry might be of special physical interest. With regard to this, 
VA leads to a simple expression for <x m ; n (we just expand Eq. ([55]) in a Taylor series around J = 0): 

CTmin ~ 2VJ . (36) 

The essential feature in Eq. (B6[) is a square root dependence. To comply with the strict analytical result for the 
unperturbed case, see Eqs. (^3|) and (E?]), we adopt 



For N = 3, a similar dependence is found to be 



VJ, (37) 



« yJyfiJ- (38) 
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as an approximation. The bold solid lines corresponding to these approximations are shown in Fig. |3J and Fig. |4|. 

Apparently, the smaller the dissipation, the better the approximations ([37]) and (|3£j| ) work. In the presence of the 
dissipation, the fluxons can exist at velocities somewhat larger than in the idealized model. The latter circumstance 
is quite natural, as in the single. Josephson junction the onset of the fluxon's instability past the Swihart velocity is 
also delayed by the dissipation!-!. 

We can also obtain an analytical approximation for N — ► oo. In this case, 

a min «V2J. (39) 

The stack with the uniform critical current distribution, therefore, has er m in = v2. The Cherenkov radiation will 
appear in the range V2 < er min < 2. 

It is convenient to present the results of our calculations as the plot -Waxd-SI)- Such a plot for N = 3 and = 2 1 / 4 
(uniform stack) is shown in Fig. ^. The region where the fluxon moves faster than c_ is shaded. It is the domain of 
existence of the Cherenkov radiation. The range of S in Fig. |^ corresponds to the maximum value of the coupling 
parameter l^jmax in 3-fold stack which is equal to l/y/2. 



V. CONCLUSION 



We have demonstrated that a single fluxon moving in the stack of Josephson junctions has the maximum velocity 
which does not necessarily coincide with one of the Swihart velocities of the system. The dependence u max (J) was 



put forward on the basis 
and (p8[). Our results 



studied numerically for N = 2 and 3. An analytical approximation for this dependence was 
of the variational approximation in the limit of small dissipation and is given by Eqs. ( [37 
show that u m ax > c_ for J < 1, in the case N = 2, and for J < y/2, in the case N = 3. This leads to Cherenkov 
radiation of plasma waves by a fluxon in the velocity range c_ < u < it max - Simulations also show that the damping 
stabilizes the fluxon motion at higher velocities, i.e. , u max (a > 0) > w m ax(a = 0)- in uniform stack with large 
N (e.g. intrinsic Josephson stacks in crystals of high-T c superconductors), exceeds c_ and, therefore, a single 

fluxon always generates Cherenkov radiation. Experiments with high-T c stacksli3 show that flux-flow branches do not 
have vanishing differential resistance Rd close to the top of the flux-flow step as in the single LJJ but bend towards 
higher voltages, in agreement with our results. From the theoretical point of view, the absence of the "relativistic" 
singularity at u — c_ is a result of the lack of the Lorentz invariance in the coupled sine-Gordon equations. 
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FIGURE CAPTIONS 



FIG. 1. The IVC's u(-y) of a 2-fold stack with a single fluxon trapped in one LJJ (the state [ljO]) for three different values 
of J = 0.5, 1, and 2. The branches for J = 0.5 and J = 1 terminate at a value M max at which the fluxon becomes unstable 
and the system switches to high voltages. The inset is a blowup of the rectangular area on the main diagram. The simulations 
were performed for a = 0.1 and demonstrate the dependence of u max on J. 

FIG. 2. The IVC's of the 3-fold stack with a fluxon in the central layer for three different values of a = 0.02, 0.04, and 0.1. 
The simulations performed for J = 1 demonstrate the dependence of it max on a. 

FIG. 3. The dependence of cr^ on J for three different values of a = 0.02 (circles), 0.04 (squares), and 0.10 (diamonds). 
The bold solid line shows the analytical dependence ( pS?] ) produced by the variational approximation. 
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The dependence of on J for three different values of a = 0.02 (circles), 0.04 (squares), and 0.10 (diamonds). 



FIG. 4 

The bold solid line shows the analytical prediction 



(p8j) based on the variational approximation. 



FIG. 5. The dependence of maximum velocity M max and Swihart velocities c± on the coupling strength \S\ for three coupled 
junctions at J = 1. The shaded area shows the domain of the Cherenkov radiation. 
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TABLE I. The values of cr m in for different numbers N of junctions in the stack. As N — > 



8 




velocity, u [norm, units] 



E. Goldobin etal. Fig. 1 




velocity, u [norm, units] 



E. Goldobin etal. Fig. 2 




ratio of critical currents, J 



E. Goldobin etal., Fig. 3 




E. Goldobin et al. Fig. 4 




E.Goldobin etal. Fig. 5 



